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Abstract 

The method of efficient description of long-term behavior of solutions of 
the non-linear first order ODE 0 + sin 0 = / for arbitrary periodic / is 
discussed. The criterion enabling one to separate and identify the quali¬ 
tatively different solutions is established. The applications of the method 
to the modeling of dynamics of overdamped Josephson junctions in super¬ 
conductors are outlined. 


Preliminaries 

A simple model of a Josephson junction [I] proposed by Stewart and McCumber 
01 is based, from the mathematical point of view, on the non-linear second 
order ODE 

00 + 0 + sin 0 = /. (1) 

Here 0 = 0(f) is the unknown function (called phase ) representing the difference 
of phases of the so called order parameter which plays the role of the macroscopic 
wave function describing the states of two weakly coupled superconducting elec¬ 
trodes constituting Josephson junction. The connection of the phase with ob¬ 
servable quantities is established by the Josephson equation which connects 
the voltage across the junction with the time derivative of 0: 


$0 dcf) ■ _ dcf) 
27r dr dt 


( 2 ) 


Here the physical constant <f>o = h{2e)~ x ( h is the Plank constant, e is the 
electron charge) is called the (magnetic) flux quantum, the parameter f used 
in Eq. m is the rescaled dimensionless time f = a; c r, where r is the ‘genuine’ 
(dimensional) time, the constant tv c called characteristic frequency is connected 
with junction properties. It can be defined by the equation u> c = 27r$Q 1 R N I C , 
where Rn, I c are the junction characteristics, namely, R^ is its resistance in 
the normal state, I c is called critical current (it estimates the maximal current 
which can flow through the junction without application of external voltage). The 
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positive constant (3 is connected with the junction capacitance, the less the latter, 
the less the former. The function f = fit) represents the current (normalized to 
I c and named bias) circulating in the junction circuit and supplied by an external 
source, ft is assumed to be specified in advance. 

Here we consider the case of a periodic bias of the known period T and 
arbitrary profile. Thus / is assumed to satisfy the condition 

m = f(t+T). ( 3 ) 

Usually / is a continuous (or smooth) function but a finite number of finite jumps 
on the period interval is also allowed 1 . 

ft is worth noting that there is a common practice in the physical literature 
to divide the periodic bias function / into some constant constituent (‘direct 
current’, DC) tdc and the residual one t ac (alternating, high frequency, rf current, 
AC): 

f = L dc + htc- (4) 

Following this convention, we assume for definiteness that the average value of / 
is assigned to tdc- Accordingly, by definition, the average value of the residual t ac 
vanishes 2 . 

The important property of a Josephson junction, which is perfectly captured 
by the Stewart-McCumber model in spite of its comparative homeliness, is the 
phase-locking effect. Qualitatively, the phase-locking may be regarded as the 
relaxation in the course of the phase evolution of the ‘loose constituent’ of the 
phase function reflecting its initial state. Then the only surviving dynamical 
‘process’ proves a steady one ‘dragged’ by the periodic bias and not depending 
of the past phase evolution. Said another way, once some lapse depending on the 
specific rate of relaxation of initial perturbations and the ‘magnitude’ of the latter 
elapses, the phase function profile proceeds with the reproducing itself , modulo a 
27r-aliquot contribution, on each subsequent time step of duration T. Obviously, 
such a situation would take place, in particular, if, regardless of the initial state 
the phase starts to evolve from, 0(f) converges for large t to a periodic function of 
the period T. There are also other, non-periodic (but still fairly specific as we shall 
see below) forms of such an asymptotically ‘steady’ phase evolution. On the other 
hand, the phase-locking property is in no way a universal one to be automatically 
attributed to solutions of Eq. dU . Depending on the model parameters, the phase 
evolution may be completely different, revealing, in particular, no periodicity or 
any other apparent order [S]. 

having added of a finite number of Dirac <5-likc pulses, the problem retains meaningful but 
such a generalization will not be pursued here. 

2 Such an interpretation reflects the typical experimental setup which fixes (, ac and considers 
id c as a free parameter to be varied. In particular, the current-voltage (I-V) curve of a Josephson 
junction is understood just as the dependence of the average voltage across it on the DC bias 
contribution whereas (, ac is kept unchanged. 
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Ph 1 (vertically uniformly shifted) 



Figure 1: Three subsequent segments (red, brown=l (red+green). green) of the sin¬ 
gle graph describing evolution of the phase function on the time interval [0, 3T] are 
placed over the common t -axes segment [0, T] by means of appropriate ‘horizontal’ 
shifts ‘backwards in time’ (plot segmenting). The specific parameter value is idc = 0.8 
the values of other parameters are given in the main text. In order to visually resolve 
the neighbor graphs here (and in subsequent similar plots) all the curves, except for a 
starting one, are uniformly shifted in the ‘vertical’ direction by an additional amount 
of space which is uniformly incremented for each subsequent graph curve (here by 0.04 
units). As one can see, the third (green) graph segment already coincides with the 
second one, essentially, yielding therefore the (approximate) specimen of the limiting 
phase evolution function which proves here periodic. 
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The plots shown in figure Q illustrates the simplest case of the phase-locking in 
a numerical example. Here the result of a straightforward numerical integration 
of Eq. (P) is displayed. For the sake of definiteness, the bias function f(t) is chosen 
to be the periodic rectangular pulse sequence, one of the pulses being shown in 
the inset. Its integral magnitude amounts to 1.5, the frequency is 27r/T = 0.47 
( T ~ 13.3685), (3 = 0.02, the width of the pulse peak (here being rather a 
plateau symmetric with respect to the plot center) amounts to 20% of the total 
pulse period. The ‘constant bias constituent’ l& c (the average of /) equals 0.8. 

The phase evolution displayed in Fig. |Tj starts with the null initial conditions 
0(0) = 0 = 0(0) chosen for simplicity reasons 3 . The figure displays the phase 
evolution during 3 sequential steps of variation of t, each of duration T. The 
graph of the first period phase evolution is of red color. The second period plot 
graph ( t E [T,2T]) is brown. To be placed over the same t- axes segment as the 
first graph, it is shifted to the left (‘backwards in time’) just by T. Similarly, 
the graph of the phase function over the third period time step (t E [2T, 3T], the 
green curve) is shifted to the left by 2 T that places it over the same t- axes segment 
[0, T], As a result, all the three subsequent portions of the phase evolution graph 
are displayed over the common segment of the f-axes becoming more eligible for 
a visual collation and the monitoring of the T-scale periodicity. 

In passing, it is worth emphasizing that due to periodicity of / the functions 
resulted from the above ‘translations’ by means of the f-shifts aliquot T retain 
to verify Eq. ©• Thus, applying such a trick, the ‘contiguous’ phase function 
representing the phase evolution over some time interval, perhaps semi-infinite, 
can be equally well represented by the sequence of phase functions (solutions of 
Eq. (£Q)) each defined on the finite interval [0, T\. 

Besides, for the better clarity, some additional artificial uniformly increment¬ 
ing vertical shifts are applied to separate graph segments in Fig. |T| The point is 
that we would like to recognize the limiting curve to which the sequence of the 
displayed ones converges, and the details of such a conversion, provided it takes 
place. Obviously, it is reasonable to add some uniformly incremented vertical 
space to the vertical positions of separate graph segments in order to visually 
distinguish the far mutually close curves which otherwise would be seen in the 
plot overlaid in a messy way. Specifically, in the case of Fig. [TJ the vertical shifts 
amounting to additional 0.04 units per each subsequent graph segment are ap¬ 
plied. This allows one to visually distinguish the second (brown) and the third 
(green) phase graph segments showing simultaneously that the corresponding 
functions are very close. In subsequent similar plots, the benefit of the trick will 
be ever more clear. 

It turns out that under the conditions assumed, the phase function on the 
third time step (the third period) apparently coincides with the one on the sec- 

3 It is worth noting that the condition 0(0) = 0 is not the most natural choice since it would 
lead to the ill posedness of the Cauchy problem in the limiting case of small (3. 
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(vertically uniformly shifted) 



Figure 2: The voltage across junction (in arbitrary units) is displayed for the phase 
functions shown in Fig. |2 In order to resolve the convergent series of graphs, the 
graphs except of the starting one are uniformly shifted in the ‘vertical’ direction, the 
trick similar to one applied in Fig. ^ 

ond step and thus already represents, with accuracy characteristic for the plot 
resolution, the desirable limiting (asymptotic) phase time dependence. The T- 
scale convergence proves here extremely fast. It is also obvious that in view of 
the convergence of subsequent graph segments to the common limiting curve the 
values of 0 on the left and right graph boundaries of its domain coincide, making 
evidence of the periodicity of the limiting solution of Eq. © we deal with. 

It is also of interest to consider the graphs of the corresponding voltage which 
coincides, up to the constant factor 4 , with the derivatives of the functions dis¬ 
played in Fig. QJ These are shown in Fig. [2l the vertical and horizontal shifts (the 
part of the plot segmenting trick) similar to ones in Fig. being applied. The 
voltage functions also rapidly converge to a periodic limiting one. 

It has to be noted that the phase-locking does not necessarily claim for the 
asymptotic phase function 0(f) to be periodic. Indeed, keeping in mind the 
physical aspect of the problem, the actually registered voltage across Josephson 
junction is not momentary but is efficiently averaged over many bias periods T. 
For the sake of simplicity, let us consider the result of the voltage averaging over 
a single period (the shortest time interval making sense in the situation under 
consideration), treating the averaging as the plain integrating followed by the 
division by the length of the integration interval. Then in accordance with Eq. 

4 In order to evade introduction of dimensional numerics which are irrelevant in the present 
discussion, we shall use throughout instead the function V as it is defined by Eq. m the ‘voltage 
function’ 0. The plots of 0 will be referred to below as V plots in arbitrary units. 
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m the averaging started at the moment t has to yield 

/ t+At 

V(t')dt' oc Af _1 [0(t + At) — 0(f)] (5) 

with At = T (the proportionality relation means here the dropping out of a 
dimensional universal constant). Thus the (asymptotic) periodicity of 0 means 
(V)t = 0 implying that no average voltage across the junction is observed 5 . 
However, instead of / periodicity, one may also assume that the increment of 
0 across the time step of duration T is independent on the moment t when 
the averaging starts: 0(t + T) — 0(f) = constant, asymptotically (for large but 
otherwise arbitrary t). This constant may not be arbitrary. The condition above 
is consistent with Eq. El if, asymptotically, 

0(t + T) — 0(f) = 27 rk, ( 6 ) 

where k is some integer, positive, negative, or zero. This number is named the 
phase-locking order. Evidently, the periodic phase function is the particular case 
corresponding to the zero order. Allowing arbitrary integer k, the average voltage 
proves quantized assuming, in principle, any of discrete uniformly distributed 
values 6 . 

Numerical integration of Eq. m allows one to easily confirm the feasibility of 
non-zero-order phase-locking phase functions. In particular, Fig. |3] displays the 
example of first order phase-locking. Here the only difference in the parameter 
setup with the zero order case (figures [TJ |2J) is the value of the parameter i dc 
now amounting to 1.1. The limiting phase function is here not periodic but its 
increments on the time steps of duration T apparently tends to 2n. On the other 
hand, subtracting the linearly growing contribution which insures this secular 
phase incrementing, a periodic function is leaved. 

Similarly, the further increasing of t dc yields the second order phase-locking 
phase evolution. It is displayed in Fig. 0 Here tdc = 1-4, the other parame¬ 
ters being kept unchanged. Increasing t dc , the larger order phase-locking phase 
functions can be produced as well. 

It is seen that in the cases of higher order phase-locking, the convergence of 
the sequence of segmented phase functions to the common limiting function is 
slower. Indeed, in the case of zero order phase-locking (Fig. [TJ t dc = 0.8) the third 
iteration apparently coincided with the second one and yield, in fact, the limiting 
curve (the phase-locking order can be estimated from the asymptotic relation 

k ~ (27r)” 1 [0(f + T) -0(f)] (7) 

5 One of the forms of the Josephson effect: some [averaged] bias current flows through the 
junction but the [averaged] voltage across it is zero. 

6 More generally, the average voltage quantization may also occur if k is a rational number, 
provided the averaging is carried out over sufficiently large number of periods. 
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Phi(vertically uniformly shifted) 



V (vertically uniformly shifted) 



Figure 3: Example of the segmented phase function (top panel) and segmented voltage 
function (lower panel, arbitrary units) for the first order phase-locking state. The 
parameters of the model are the same as in figures |T| |21 except for id c = 1.1. The 
color of subsequent phase function segments uniformly varies from red to green. The 
uniformly incremented shifts in the ‘vertical’ direction are applied including the first 
order phase-locking responsible contribution —2n for the phase (—27T+0.2 units for the 
top panel and 0.2 units for the lower one per subsequent graph segment, respectively). 
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(vertically uniformly shifted) 



V (vertically uniformly shifted) 



Figure 4: Example of the phase-locking of the order two. The parameters of the model 
are the same as in figures mm except for z.d c = 1.4. The color of subsequent phase 
function segments (top panel) and voltage function segments (lower panel, arbitrary 
units) uniformly varies from red to green. The uniformly incremented shifts in the 
‘vertical’ direction are similar to ones applied in Fig. 0 differing in the shift contribution 
for the phase connected with the non-zero order of the phase-locking which amounts 
here to — 47r instead of —2n applied therein. 























Phi(vertically uniformly shifted) 



V (vertically uniformly shifted) 



Figure 5: Example of asymptotically non-convergent (‘irregular’) evolution. The pa¬ 
rameters of the model are the same as in figures^Gl except of t dc = 1-45. The uniformly 
incremented shifts in the ‘vertical’ direction (— 47T + 0.2 units for the phase and 0.2 units 
for the voltage functions per subsequent graph segment, respectively) are applied. 


taking place for sufficiently large t (formally, in the limit t —> oo), provided the 
phase-locking is developed). In the case of the phase-locking of the first order 
observed in the case l ( \ c = 1.1 the four T-lapses are needed for the reaching 
of a comparable closeness to the limit. For the phase-locking of the order two 
which takes place, in particular, if t ( \ c = 1.4 the number of necessary iterations 
increases up to 15. And it suffices to set up 6 dc = 1.45 to lose the apparent 
convergence to any steady asymptotic state at all. The corresponding apparently 
non-convergent segmented phase and voltage function plots are shown in Fig. 
El These reveal no tendency to converge to any limiting curves demonstrating 
‘irregular’ (or even apparently ‘chaotic’, as far as one concerns the variations on 
the scale T) behavior. The value of the expression (171) . irregularly oscillating 
around 2, does not reveal apparent tending to any limit as well. 

It is worth mentioning that the convergence — or the divergence — of the 
sequence of phase functions describing a single long-term phase evolution is often 
difficult to reliably establish by means of straightforward numerical integrating of 
Eq. (El) . Indeed, the slower such a sequence converges, the longer fragment of the 
phase evolution has to be analyzed. Accordingly, approaching the hypothetical 
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boundary of the area in the space of the bias functions where the evolution of the 
phase leads to a steady state ( phase-locking area in the parameter space), one 
inevitably encounters with insufficiency of the available computational resources. 
The very boundary, where the convergence reverses to the divergence, is not 
detectable in this way. Besides, in view of the obscurity of the corresponding 
time scale, the ‘strong’ phase divergence (more exactly, the absence of convergence 
to any asymptotic steady state) is even harder to detect through the apparent 
behavior of numerically generated phase function than the too slow convergence. 

The formulation and substantiation of the method allowing one to efficiently 
distinguish the phase-locking property in the important particular case of ‘over¬ 
damped’ version of equation m is the main theme of the present discussion. The 
term ‘overdamped’ refers here the case 

0 < 1. (8) 

The parameter scales the (only) term involving the second order derivative of 
the phase function in Eq. ©• It might be supposed that under certain conditions 
including imposing of limitation on the second derivative values the relevant 
solutions of Eq. © may be approximated by the solutions of the equation 

0 + sin 0 = / (9) 

obtained from Eq. © by a mere discarding of the second order derivative term. 

The numerical example illustrating the impact of the above ‘simplification’ is 
shown in Fig. © the discussion of the status of the corresponding approximation 
for ‘overdamped’ systems (in the case of sinusoidal bias /) can be found in [5]. 

In Fig. ©the black curve represents the phase function computed in accordance 
with Eq. ©. the red one is the corresponding solution of the general RSJ-equation 
©• The ‘red’ solution is shifted downward by 0.1 units in order to allow easier 
visual distinguishing of the close graphs on the common plot panel. The green 
curve represents the difference of the two phase functions above which is magnified 
by the factor of 10. The bias pulse is shown in the inset. Its form is similar to 
the ones used above, see Fig. [TJ except for the bias ‘DC constituent’, td c , which is 
here chosen vanishing (i.e. f Q f(t)dt = 0), and the integral magnitude of a single 
pulse which here and in all the numerical examples below amounts to 3.5. 

It is seen that under the conditions assumed the solutions of Eqs. © and 
© are very close, in total, both in shape and in magnitude, the deviations 
emerging at several moments of time ultimately relaxing. Specifically, the solution 
difference reveals the splashes originated in three points: at t = 0, and on the 
front and back ‘shocks’ of the bias pulse. At right near them, the difference 
soon reaches a maximum and then quickly die. In total, the distinction of the 
corresponding solutions of Eqs. © and © may be considered fairly mild. 

It is easy to see that the splash following the point t = 0 arises due to the 
specific choice of the initial conditions for solution of Eq. © which read 0(0) = 


10 



Figure 6: The related solutions of Eqs. 0 (the red graph) and to 0 (the black one) 
for the same periodic rectangular pulse bias (shown in inset) displayed together with 
their difference magnified by the factor of 10 (the green graph) . The red plot is shifted 
downward by 0.1 units. The bias is the same as above except the parameter c which 
is now assumed to vanish and the pulse magnitude here equal to 3.5. 


0 = 0(0). Then Eq. 0) implies 0(0) = (5 1 /(0) which, for /(0) ^ 0, is a large 
quantity contrary to assumption required for approximating Eq. 0 by 0. In 
particular, it leads to the fast growth of the first derivative immediately after the 
start of the phase evolution. This effect is also observed in Fig. [21 

Apart of the apprehensible effect of non-equivalent initial conditions, the cause 
of the other discrepancy splashes is also quite understandable: they arise due to 
the discontinuities of the R.H.S. function /. Indeed, as a pulse front/back jump 
occurs, the discontinuity in / converts to a finite jump in the derivative 0 in 
solution of Eq. © which, in turn, emerges the Dirac 5 function-like irregularity 
in the second derivative 0 which again violate the assumption of limitation of 0 
magnitude required for the approximating of solutions of © by solutions of Eq. 

©. 

Evidently, the cause of these peculiarities, leading to major distinctions of 
solutions of Eq. 0 and Eq. 0 shown in Fig. [HI originates in either the admitting 
somewhat inadequate initial conditions or in a too rough representation (by a 
discontinuous function) of the bias term. Since for a more physically realistic bias 
model no h-like irregularities of the phase derivative can arise, the agreement of 
the models based on Eqs. 0 and © would improve. 

As compared to more deep Eq. 0 whose study is still possible by numerical 
methods alone, Eq. 0 proves more transparent and allows fairly deep analytical 
treatment. It seems also to be of interest in its own rights as the example of a 
non-linear problem of a deep physical relevance associated with a simple linear 
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dynamical system 7 . Eq. © possesses therefore a number of remarkable properties 
which are the subject of the discourse below. 


The modeling of the phase dynamics 

Master identity and its consequences 

The approach to the problem of description of generic properties of solutions 
of Eq. @ on the timescale exceeding the period T including their asymptotic 
behavior can be based on the following observation: 

Master identity 1 For any piecewise smooth continuous functions <f> = (j>{t), 0o = 
0o (t), let us define the functions 


Po 

= Po(t) = P o [0o](i) = J COS0O (t)dt, 

(10) 

Qo 

= Qo{t) = Qo[0o](^) = j e~ Po{t) sin0 o (t) dt, 

(11) 

and introduce the notations: 


c 

= C(*) = e 1 *<‘>, Co = Co(() = 

(12) 

C 

— C[0> <t>o\(t) — Qo + ie , 

(13) 


Then the following identity takes place 

1 dF ~ 1 

-K-Co^^^Co-DKl-C-DKo] (14) 

with the differential operator D[-] is defined as follows 

D IA = + ( 15 ) 

for some function f = fit ) and arbitrary piecewise smooth z = z(t). 

Its proof reduces to a straightforward computation applying definitions of the 
functions involved.□ 

Remarks: 

' In the case of the bias function f{t) = B + Acosuj{t — t 0 ), where A,B,u,to are some 
constant parameters, which is most important from viewpoint of applications, Eq. © proves 
converting to the following remarkable equation: 

{« 2 - S 104 + B)c!+ - 4 - ff ]^ + ^ 2 } T = 0 
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• For the sake of definiteness, the indefinite integrals in Eqs. (fTUITTD can 
be understood as ... d t. Another choice of lower integration boundary 
would lead to a linear transformation of C\.. .] _1 : an overall constant factor 
is induced by the changing the lower boundary in the P-integral while a 
change of the lower boundary in Q-integral definition yields a constant 
summand. 


• Solving Eq. m with respect to £, one gets the equation 


l + C(Q 0 -ie~ p °) 
1 + C(Qo + ie _p °) 


(16) 


Assuming P 0l Q 0 to be real, for complex valued Co, C both of these variables 
can be unimodular, |Co| = 1 = |C| (and (po,(j) real, see (TH21 A if and only if 
C is real. 


The connection of the relationships inferred from the identity (ED with the 
equation © follows from the identity 

^ + ^(C 2 -l)-i/C = iC(0 + sin0-/), (17) 

provided ( and are connected by Eq. ED- In conjunction with master identity, 
it immediately yields the following j7J 

Proposition 1 Let <f> Q (t ) verify Eq. (0J) whereas P 0 (t),Q 0 (t ) be defined by Eqs. 
W, HI- Then any solution <f>[t ) of Eq. (OJ) satisfies either the equation 


e m) _ c iMt ) 1 + C (Qo- [e P °) 

1 + C(Qo + ie~ p o) ’ 

where C is some real constant, or the equation 


gi <!>{*) — e i0o (*) 


Qo ~ ie p ° 


Q 0 + ie p o ’ 

the latter being in fact the limiting form of the former as C 


(18) 


(19) 


oo. 


Remarks: 


• When regarded as the relationship determining 0, Eq. ED or ED spec¬ 
ifies it not uniquely but modulo 2n. On the other hand, this is the only 
arbitrariness involved in such a definition of (fit). 

• The initial, ‘ground’ solution 0 O is produced by the same equation C3 in 
the case C = 0 when ED converts to the identical transformation (modulo 
2n). That is why just C rather than C~ l which has a formally simpler 
representation in terms of 0’s, see ED, is employed in ED- 
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• Given 0(0), the constant C can be easily found. Indeed, it follows from Eq. 
m that, placing the left boundary of the integration interval at t — 0, one 
has -Po(O) = 0 = Qo(0), that yields 

A<H0) _ p i<M0) 

r = c = 

gi^(O) _|_ gi</>o(0) 

Thus, given any particular solution of Eq. © and the integrals P 0 (f)> Qo(t), 
defined from it by Eqs. m, d and obeying the conditions P 0 (0) — 
0 = Q 0 (0); the solution of the Cauchy problem for Eq. © is explicitly 
represented (modulo 2n) by Eqs. (fl8l) and da 8 . 

• In view of PI) one also gets 

C[0Oj0]|t=O = —G[0, 0o]|t=O (21) 

Fig. 0 illustrates the above statements. Here the top panels contain the plots 
of the two different solutions of Eq. ©. The bias function / is the same sequence 
of periodically repeated rectangular pulses which was used above, see the inset 
in Fig. El the only difference is the constant constituent which now amounts to 
6dc = 1-45 (it is worth reminding also that now (3 — 0). The phase function of 
the first period phase evolution, starting at the coordinate origin, is considered as 
the ‘ground’ solution of Eq. © denoted above as 0 o- If is plotted at the top-left 
panel. In the lower-left panel, the integrals Pq, Q 0 calculated from 0 O are plotted. 
The top-right panel displays the ‘continuation’ of 0 O to the second period [T, 2 T] 
which, afterwards, is returned (‘shifted’ to the left by T) to the 0 0 domain. (No 
‘vertical’ shift is here applied.) As it had been mentioned above, due to the 
periodicity of /(f), the ‘shifted’ solution is again a solution of Eq. ©. We denote 
it 0(f). (This function is distinguished by the specific initial condition 0(0) = 
0o(T) which however plays no role in the current context.) Finally, the graph 
displayed in the lower-right panel is the result of straightforward computation of 
C = C[0, 0 o ](f) in accordance with definition (THU) and <H2D- More exactly, the 
deviation of C off its averaged (on the interval [0, T]) value which is Ri—0.406445 
is shown. 

One sees that the functional C = G[0,0o](f) computed ‘from the first prin¬ 
ciples’ is, after all, a constant up to a small deviation. The latter irregularly 
oscillates around zero with the amplitude which is at least 6 orders less than the 
phase magnitudes. It represents the ‘numerical noise’ reflecting mostly tolerable 
inaccuracy of approximate numerical solutions of differential equation. 


[0(0) - 0 o (O)]. (20) 


8 Although the case 0(0) — 0o(O) = 7 r mod 27 t is not covered, formally, by the above formulae, 
a minor obvious modification corresponding to the transition to the limit C —> 00 allows to 
treat it as well. 
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Figure 7: The top-left and top-right panels display the phase functions 4>o(t), 
respectively, verifying Eq. ©• At the lower-left panel, the integrals Po{t) (the red 
curve) and Qo(t) (the green curve) are plotted, see Eqs. (1101) . (ITT1) . The deviation 
of C[(j), 4>o\(t.) computed in accordance with Eqs. (THU) . (TH21) from its average value is 
displayed at the lower-right panel. 
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(with some vertical shifts) 



Figure 8: The black curve represents the ‘ground’ phase function. The red and green 
curves are the phase functions (shifted downward by 0.1 and 0.2 units, respectively) 
which are computed in different ways: the red curve is obtained by means of Eq. m 
for C = 0.2 while the green one is the result of numerical integration of Eq. m with 
the same initial conditions as the ‘red’ solution obeys. Finally, the blue curve shows 
the difference of the “red” and “green” functions magnified by the factor of 0.5 x 10 6 . 

The closely related computation displayed in Fig. |H1 illustrates the application 
of Eq. (1T%1) for the generating of new solutions of Eq. © from a known one. Here 
the black curve shows the ‘ground’ solution, 0o, starting in the coordinate origin; 
note that any other phase function might be used instead. Then we chose, say, 
C = 0.2 and compute <f(t) by means of Eq. (JTHJ) - The result is represented by the 
red curve which was shifted downward by 0.1 units for the convenience of further 
visual collations. Next, we carry out the straightforward numerical integrating 
of Eq. © adopting the same initial conditions <f)\t=o = 0(0) as the just generated 
solution obeys. The green curve represents the result of the integrating which is 
also additionally shifted downward , here by 0.2 units. Finally, the blue graph is 
the difference of the results of computation of the same phase function by these 
two methods (viz the application of Eq. m and the numerical ODE integrating) 
magnified by the factor of 0.5 x 10 6 . Again, one sees that Eq. (fTHl) ensures the 
stable accuracy about six true decimal digits with the discrepancy falling in the 
level of ‘numerical noise’. 

ft is important to note that, in the relationships above, the roles of the func¬ 
tions 0o and 0 (any two solutions of Eq. ©) should be symmetric by general 
reasons. However Eq. (CHI) does not reveal, apparently, such a symmetry since it 
involves the integrals Pq, Q q determined by the solution 0 O but no similar con¬ 
tribution connected with 0 is present. Besides, Eq. m is easily solvable with 
respect to 0 but represents to a nonlinear integral equation with respect to 0 q. 

Any inconsistency does not arises here however and the asymmetry mentioned 
above is actually a fallacious one. The point is that there is a specific algebraic 


16 












relation which constrains P- and Q-functionals associated with any two solutions 
of Eq. It enables one, in particular, to represent any P,Q as a simple ele¬ 
mentary function of another P, Q-pair. Namely, the following remarkable relation 
takes place. 


Proposition 2 Let(ft 0 (t), <ft{t), P 0 (i), Qo(t ) be as in Proposition QJ. Let also P(t), 
Q(t ) be determined from <ft(t) in the same way as P 0 (t), Qo(t) are determined from 
(fto(t). Then 


p _ (Qo + ie- f Q) - C 
V 1 + C(Q a + ie-«>) 


( 22 ) 


with the same real constant C. 


Remarks: 


• Inverting ( 1221 ) . one gets the alike formula 


Qo + ie P ° 


(Q + ie" p ) + C 
1 -C(Q + ie~ p ) 


(23) 


which differs from (1221) . apart of interchanged roles of P, Q and Po,Qo, by 
the opposite sign of the (7-constant alone. 

• From viewpoint of the induced transformation group structure, the trans¬ 
formation © coincides with the relativistic velocity addition rule. 


• In view of the relationship above, it seems natural to incorporate the real- 
valued functions P(t), Q(t) into the complex-valued one defined as follows: 

P(t) = Q(t) + ie- p(t) . (24) 


Then the transformation described Eq. (1221) may be named, in a sense, 
meromorphic since P is expressed via P 0 as a (meromorphic) function of 
Po- 

• The definition of P, Q-integrals by m, © -like equations is equivalent to 
the equation 

Ap(£) = _i eW^pm. (25) 

at 

Indeed, rewriting it as d(3?P)/ dt + id(3P)/ d t — (sin (ft — i cos (ft) (OP) and 
separating the pure imaginary part, one gets d(QP)/d t — —cos (ft ■ OP, 
which, for nonzero OP, yields OP = exp (— f cosc/>d</>), i.e. in accordance 
with definition (1221) . Eq. m in fact. Next, separating the real part, one 
gets d(9tp)/ dt — sin (ft ■ OP. Integrating it and taking into account the 
representation of OP just obtained, one gets Eq. (ITTlt .D 
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• Working in terms of P instead of P, Q, the choice of t = 0 as the lower 
integral limit in (HOD, (ED, is equivalent to the initial condition 


T(0) = i 


(26) 


for the function P verifying Eq. 

• Inverting Eq. (1221) . P 0 can be represented as a function of P (see (1231) 1. 
Then Eq. © can be solved with respect to 0 O (t) which is represented as 
explicit function of 4>(t), P(t), Q(t ) (and C) as follows 


e iMt) = C W ) 1 Z Z ie ) 

' 1 — (^(Q + ie _p ) ’ 


(27) 


The dual Eqs. ED,USD manifest the symmetric role of the two solutions 
d>, 6 n noted above. 


Proposition proof. The relationship to be proven belongs to the category of ones 
for which, as the true formula is recorded, the proof reduces to a straightforward 
computation. Indeed, let us rewrite Eq. eg as follows 


e i 4>C) 


= i0o(t) 1 + CPo 
1 + CP o 


and calculate the t-derivative of the difference 


Po ~C r 1 1 + C 2 

i + cp 0 c c(i + cp 0 y 


Then the straightforward computation yields 


d<5 
d t 


dP l + C 2 dP 0 
~dt (1 + CP 0 ) 2 d t 


2 C 


i<t> 


(P - P) - 


1 +c 2 

(1 + CPof 


d<t> o 


(P o — Pq 


i p i0o 

2 C 


fl + CPo 
\1 + CP 0 


(P-P) 


(1 + CPo) 2 ° 



The following identity takes place 


P-P 


1 + C 2 f l 1 \ 

C \1 + CPo ~ 1 + CPo) 


5-5 + 


1 + C 2 


(1 + CPq)(1 + CPq) 


=— (P 0 - P 0 


(28) 
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It implies, together with the equation just derived, the following series of equali¬ 
ties: 


d 6 
d t 


2 


2 


2 


ig^o 


ig^o 


ig^o 



The last equation is equivalent to the system of two linear homogeneous first 
order ODEs for the two (real valued) functions JM, 5s5. Moreover, in accordance 
with standard ‘normalization ’ m and 5 definition, one has 


<5(0) = 0, 


(29) 


the null initial condition. Thus S(t) = 0, i.e. one gets the equation 


which is nothing but Eq. (1221) . □ 

Thus, in accordance with the remarks above, we may consider the expression 
m as the functional over the set of pairs of solutions of Eq. @ (the direct 
square of the solution set) taking it onto the real axis. To be more exact, since 
the limiting case C —> oo is quite legitimate, one has to replenish the real axis by 
the infinitely remote point. The result is the circumference but it appears here 
as the isomorphic image of the real projective line RP 1 which is just the natural 
‘index space’ to be used for the ‘enumerating’ of solutions of Eq. ©• Specifically, 
having fixed some 0 O , the inverse map from S 1 ~ RP 1 onto the Erst factor in 
the direct product yields the (1-1 non-canonical) parameterization of the space 
of solutions by the circumference points. 

It is worth noting that the (7-map m, m is antisymmetric with respect to 
the interchange of the functional arguments (j), 4>o, he. 


C[(t>A o] 


C[(j) 0 ,(j)\. 


(31) 


The property m can be established as follows. Eq m implies 
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and therefore C[(f> o, 0] is automatically a constant, provided C[0, 0o] is- Fur¬ 
ther, in accordance with m the equality —C[0o,0] — C*[0, 0 O ] takes place for 
t = 0. Hence these constants coincide up to the opposite signs and CD holds 
everywhere. 

Proposition 3 If 0, 0o are solutions of Eq (© then Eq. m is satisfied. □ 

Phase-locking and its necessary condition 

Now let us consider the implications of the above relationships in application 
to the property of the phase-locking. The latter is connected with the specific 
asymptotic behavior of the phase functions 0(t) verifying, in our case, Eq. ©. 
Namely, in the case of phase-locking Eq. © has to be satisfied, asymptotically. 

To describe efficiently this property, let us define the sequence of functions 
<j>j(t) defined on the segment [0, T] as follows 9 

0j(f) = 4>{t + jT) -27t[[0(jT)/27t]], j = 0,1,2,3..., t e [0,T] (33) 

where the double brackets [[...]] stands for the integer part of the real number 
enclosed. Thus the plot of (j>j(t) displays how the ‘genuine’ phase function 0(f) 
looks like ‘on the j’th segment’ of the length (duration) T with respect to the 
closest level aliquot to 2i r. To that end, its graph is shifted downward or upward 
by such a number of ‘full phase revolutions’ 2ir which returns its left-boundary 
point 0j(O) to the segment [0, 2tt). 

We have already used such a trick in the arrangement of plots of phase func¬ 
tions calling it occasionally ‘the segmenting’. Here the explicit transformation of 
the phase function on the corresponding t-segments of duration T yielding a spe¬ 
cific sequence of the phase functions defined on the segment [0,T] is introduced. 

All the functions <f>j satisfy Eq. ©- The specific (and characteristic) property 
of the sequence of its solutions 00s associated with a single solution 0 defined 
for all t is obviously the following: 

(f>j+ 1(0) = 0j(T) mod 271,0^(0) G [0,2vr),j = 0,1,2... (34) 

Given the sequence {0j} of solutions of Eq. © on the segment [0, T] fulfilling the 
condition (ETA . the valid phase function 0(t),f G [0, oo), can be reconstructed in 
the obvious way. 

9 There is some abuse in these notations prone of a mix of the meanings of the symbol 0o- 
Above, it was understood as one of the phase functions from the pair 0,0o (see e.g. JT5J) or as 
the ‘ground’ solution obeying the condition 0 o (O) = 0 as in Fig. [3 such a usage leading to no 
interpretation problems. The new interpretation of the same symbol is the particular element 
of the sequence of functions defined by Eq. (pf . It does not match the above. Usually, it is 
clear what is meant. However, if below the both meaning loads of the symbol 0o (as well as 
Po, Q o) meet in a common context, a slightly modified notation, 0(o), will be employed for its 
first interpretation. 
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It is convenient to represent the phase-locking property in terms of properties 
of the functional sequence 0 r It can be stated the following: 

In the case of phase-locking the sequence (f>j uniformly converges on 
the segment [0, T] to some limiting function 

0oo (t) = hm 4>j(t) (35) 

j-»oo 

which satisfies Eq. © and the condition 

(T) = 0^(0) + 2vr k e^ (T) = e^ (0) (36) 


for some integer k. 


It is natural to select and fix the ‘ground’ specimen of phase functions denoting 
it 0 o (t) , t e [0,T], characterizing it by means of the null initial value 0 o (O) = 0 
(the symbol 0( O )(f) and its derivates will be also used in cases prone of confusion 
with the notation introduced in (|S3D). Having posed the problem in such a way, 
0o (t) is completely determined by the bias function /(f). We also assume the 
value t = 0 to be the lower integration boundary for the integrals in the P 0 , Q 0 
definitions (fTO . (fTTl) . Then Eq. (fTHl) implies that there exists a constant such 
that 

c i<foo(t) _ c i0o(*) 1 + Coo(Qo(^) ~ le P ° W ) (gyN 

l + CM^+ie-roUy ^ 

(We assume to be hnite but the case of ‘infinite C^ is also tractable, with 
minor modiheations, in the same way.) Taking the above representation of e 1 ^ 00 ^ 
into account, it is easy to see that Eq. (dD is satisfied if and only if the equation 


cl x (Qo(T) cos i MT) + e- p » (T) sin jfo(T)) 

+c a o x (Qo(T) sin l<po(T) + (1 - e~ mT> ) cosJifo(T)) 
+ sin i<p 0 (T) = 0 


(38) 


admits a real root C^. In turn, the latter condition is equivalent to the inequality 
A[/] = -Q 0 (T)(e“ Po(T) + l)sin0 o (T) 

-I (-Qo(T) 2 + {e~ Po ^ + l) 2 ) (1 - cos0o(T)) 


+ (1 


_ p -Po(T)\2 


4e _Po(T) x 

_i e h p o( T )Q 0 (T) sin |0o(T) + cosh i P 0 (T) cos |0 O (T) 


- 1 


> 0 


(39) 


(40) 




cosh |P 0 (T) cos §0o(T) - ie 2 Po(T) Q 0 (r) sin §0 O (T) 


> 1. 


(41) 
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The above is therefore the necessary condition of phase-locking. It reflects in fact 
the property of the very function /. 

As we shall show, the similar but strict inequality is the sufficient condition 
of the phase-locking for the phase function described by Eq. Q. (In the case of 
the equality, the convergence to some asymptotic limit is observed as well but 
it is slower and reveals some other specialities.) The vantage of such a form of 
the phase-locking ‘monitoring’ is that this is the asymptotic property manifesting 
itself for sufficiently large t. Generally speaking, one cannot forecast in advance 
how long phase evolution has to be tracked for the detection of the corresponding 
T-scale reproducibility of the phase function form which would make evidence of 
the phase-locking. On the other hand, making use of the condition (14111 . the 
appearance of phase-locking is deduced directly from the properties of a single 
solution of Eq. @ computed on the finite interval [0, T\. 

More on the role of the discriminant A[/] 

The function A = A[/] plays an important role in the problem of description 
of asymptotic properties of solutions of Eq. (JH- We may interpret it as the 
functional on the set of bias functions / since the functions (poi Po, Qo from which 
it is built upon are uniquely determined by /: 0o is the solution of Eq. © with 
initial condition 0o(O) = 0, Pq, Qo are calculated from <fo in accordance with Eqs. 
EH, EH) and the initial conditions -Po(O) = 0,Qo(0) = 0. One may also choose 
/ to belong to some more restricted class of functions, for example, a family 
parameterized by a finite set of parameters. Then, one may also regard A as the 
function of these parameters and study its properties. 

Adopting the last interpretation, Fig. El shows A as the function of a single 
variable 10 , the constant bias constituent tdc, assuming the very bias function to 
be the periodic sequence of rectangular pulses shown in the inset in Fig. |H 

One recognizes the three minima on the fragment of the A plot seated in Fig. 
El the left one situating very close to the horizontal coordinate axis and the middle 
one being fairly steep. More narrowly, the minima above are shown in Fig. E3 
where the plots of their vicinities are displayed with higher resolutions. One sees 
that in all the three cases the minima lay well below the horizontal coordinate 
axes, ffence each of them is encompassed by some tdc segment where A assumes 
negative values (and this is the universal property: there are no positive minima 
of A). These are separated by segments where A is positive. 

As we shall show, each segment of positive A determines the area of the 
values of the parameter tdc giving rise to the phase-locking phase evolution of 
a common order. Such segments are directly associated with so called Shapiro 
steps — strictly constant voltage segments observed on the Josephson junction 

10 We omit here and below the ‘functional argument’ [/] of A, provided this cannot lead to a 
misunderstanding. 
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Delta[f] 



Figure 9: The dependence of the discriminant A on the constant constitution of 
the bias with the other parameters held fixed. 





Figure 10: Separate representations of vicinities of the minima observed on the graph 
shown in Fig. El are given with higher plot resolutions. 


23 
























I-V curve PHI In particular, the domain of idc values shown in Fig. El where 
A > 0, contain the steps of the orders, from left to right, —1 (extending to the 
left beyond the plot boundary), 0, 1, and 2. 

The values of i ( \ c falling into the segments of negative A correspond to the ap¬ 
parently ‘irregular’ behavior of the phase which is similar to one displayed in Fig. 
El We shall see however that these phase evolutions do not correspond to an ac¬ 
tual chaos (pseudo-chaos) [5j. Rather, these are the manifestations of a ‘beating’ 
produced by two inconsonant frequencies. In particular, such phase evolutions 
imply quite de fini te average voltages which are obtained by the averaging of © 
over large time intervals. This point will be considered in more details later on. 

Unstable and stable phase-locking candidates 

The algorithm allowing one to reconstruct the phase functions at any moment of 
time t divides into several steps. At first, one has to solve the quadratic equation 
(TTH1) with respect to C, assuming its roots to be real (that takes place iff A > 0). 
Then, if A > 0, Eq. (fTHl) yields the two phase functions which are the candidates 
to the role of possessor of the ‘refined’ property of phase-locking (the limiting 
function asymptotically approximating generic solutions). 

A somewhat delicate point is the revealing which of the two solutions is the 
one we search for. The answer is connected with the stability property (the 
phase-locking must be stable) and can be provided, in principle, with the help of 
the standard perturbational analysis (the local stability problem, cf. 0) which 
however yield not the actual problem solution but rather a recipe of its compu¬ 
tation. Fortunately, there is an attractive opportunity to directly manifest not 
only the effect of infinitesimal perturbations but, at one blow, to establish the 
‘global’ stability property allowing perturbations of arbitrary magnitudes. Here 
we seize upon it. 

This point is illustrated by Fig. E| Here the black curve shows the ‘ground’ 
solution 0o (starting from the coordinate origin) for the rectangular pulse bias of 
the profile shown in the inset in Fig. El and qic = 1-25 (which yields a positive A 
value, see Fig. ©)• 

The red curves are the graphs of solutions (‘phase-locking candidates’) one of 
which is expected to realize the ‘refined’, precisely steady phase-locking evolution 
obeying the property © and representing, in a sense, the asymptotic limit of a 
generic phase function. The phase-locking candidates are constructed by means 
of Eqs. (fTHl) with the two values of the constant C obtained from Eq. (fTRli . 

The green curves represent the phase functions constructed by means of a 
straightforward numeric integrations of Eq. © obeying the same initial condi¬ 
tions as the ‘red’ solutions does. The plot graphs are afterwards shifted ‘in vertical 
direction’ by 0.5 units upward (the green top curve) and downward (the lower 
one) making easier the visual collation of the profiles of the functions involved. 
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Phases (with some vertical shifts) 



Figure 11: The black curve shows the ‘ground’ solution </>o (starting from the coordi¬ 
nate origin). The two red curves are the the candidates to the role of the asymptotic 
limiting phase function which are constructed by means of Eqs. pg|). CHI. The green 
curves are the phase functions constructed by means of a straightforward numeric in¬ 
tegration of m obeying the same initial conditions as the ‘red’ solutions; for better 
clarity, they are shifted, afterwards, in ‘vertical’ direction by 0.5 units upward (the top 
curve) and downward (the lower curve), the red graphs being undergone no shifts. 



Figure 12: The plots display the discrepancies between the phase-locking candidates 
displayed in Fief, fill and obtained by means of Eqs. (I38[) . (fTHl) against the ones obtained 
by means of straightforward numerical integration of Eq. © with the same initial 
conditions as the former functions obey. 
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Figure 13: The evolution of perturbations of the two candidates the ‘refined’ phase¬ 
locking states represented in Fig. El by red graphs are displayed, the perturbation 
function corresponding to the stable evolution (the red graph) being magnified by the 
factor of 10 2 . 

The apparent qualitative coincidence of the results of the two ways of com¬ 
putation of ‘phase-locking candidates’ is numerically confirmed in Fig. H21 where 
their relative discrepancies are plotted. One sees that the two computations 
agree at the level one part in 10 6 which is close to the accuracy of the numerical 
integrating of the ODE involved. 

It is natural to suppose that only one of the two solutions shown as the red 
(or, equivalently, green) curves in Fig. (fTTJ) is of interest as a model of physical 
relevance because another one is necessarily unstable (it plays the own specific 
role as the separator of two ‘basins’ of the phase-locking in the space of all phase 
functions, though). Indeed, the distinction of their stability properties is clearly 
seen in Fig. El 

Here the evolution of perturbations of the two solutions shown in Fig. El 
(the red curves) is displayed, The perturbed phase functions are defined as the 
solutions 4(t) to Eq. dHJ) with initial conditions distinct from the ones the phase¬ 
locking candidates obey by amount of 0.1%: 0^(0) = 0(0) x (1 ± 0.001). For 
definiteness, we chose 4> £ ( 0) = 0(0) x (1 + 0.001) for the upper curve in Fig. W2 1 
and 0( £ )(O) = 0(0) x (1 — 0.001) for the lower one. Accordingly, the upper (green) 
graph in Fig. ITol shows 0 ( - £ '(i) — 0(t) for the ‘upper’ phase-locking candidate; it 
makes evidence of a strong instability since the deviation is permanently (and, 
in fact, exponentially) growing. On the contrary, the lower (red) curve in Fig. [HU 
demonstrate the exponential relaxation of the perturbation; moreover, to display 
it more clearly, the perturbation value displayed is magnified by the factor of 10 2 , 
i.e. the function lOO(0( £ )(f) — 0(f)) is actually plotted. 

Thus we see that only the lower red curve in Fig. El corresponds to a stable 
phase evolution. The evolution described by another, upper, (red) curve in Fig. 
El related to another solution of Eq. is unstable. As a matter of fact, 
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(with some vertical shift) 



Figure 14: P-integrals corresponding to the stable steady phase-locking are shown. 
The black curve is the P-integral for the ground phase function cj )o (the black curve 
in Fig. HU). The red curve represents P-integral is associated with the refined phase¬ 
locking candidate determined by means of Eq. (|57|). The green curve represent the same 
function obtained directly from definition m, afterwards it being shifted downward 
by 0.5 units. 


the first curve describes an attractor of all the ‘neighboring’ solutions while the 
second one ‘repels’ them. 

Above, we have numerically demonstrated the validity of the equation <mu 
allowing one to build a new solution from the known one. Similarly, Eq. m 
allows one to determine, knowing (7, the integrals P{t) and Q[t ) associated with 
the former. Figures 1771 and ITHl displays the results of the corresponding compu¬ 
tation. There the black curves represent the P- and Q-integrals corresponding 
to the ‘ground’ phase function </>o which is represented by black graph in Fig. 
HU The red curves are the integrals for the refined phase-locking candidates 
which are determined by means of Eq. CUD. The green curves (which are shifted 
in vertical direction downward) are the same functions but they are obtained by 
straightforward numerical integrating in accordance with P and Q definitions. In 
Fig. HU the relative discrepancies of the two ways of computation of the P- and 
Q- integrals are displayed. (Note that the relative discrepancy growth observed 
near the left boundary t = 0 is the numerical artifact caused by the vanishing of 
P(0) and Q(0).) 

T-scale discretization of the phase evolution 

As it has been noted, any phase function (j)(t) defined on arbitrary t domain can 
be equivalently described by the sequence of the functions 4>j{t) defined on the 
common interval [0, 27 t], see Eq. (I->MI) . At the same time, as a solution of Eq. 
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(with some vertical shift) 



Figure 15: Q-integrals corresponding to the same stable phase-locking state as in Fig. 
o are shown. The color meaning is the same as therein. The green graph is shifted 
downward by 1 unit. 




Figure 16: The relative discrepancies of the different ways of computation of P- and 
Q-integrals (by means of Eq. (1371) and the computation in accordance with definitions 
from the known phase functions) corresponding to the stable phase-locking. 


(with some vertical shift) 



Figure 17: The graphical matter displayed is similar to the one shown in Fig. [21 
but it is connected with the unstable steady phase evolution whose phase function is 
represented by the upper red curve in Fig. 




















(with some vertical shift) 



Figure 18: Similar to Fig.^Jbut corresponds to the unstable phase-locking candidate. 
The ‘ground’ (^-integral corresponding to the ‘ground’ phase function (the analogue 
to the black graph in Fig. Ej) is not shown because of the too large distinction in the 
magnitudes. 




Figure 19: The relative discrepancies of the results of two different ways of computation 
of P- and Q-integrals corresponding to the unstable phase-locking candidate. 
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Figure 20: The sequence of converging Cj which indicates the phase-locking property 
is shown. The bias is the same periodic rectangular pulse sequence as above besides 
the specific value of the DC contribution tdc = 1-46. 


©, each function <f)j(t) can be represented in the form (I Kill for the own specific 
value of the constant C which we denote Cj. A simple but important observation 
reads: 


The property of phase-locking equivalent to the existence of the limit 
4>oo(t) of the sequence of functions (p 3 (t) (see Eq. (1351) ) associated with 
a generic solution of Eq. ©, is also equivalent to the existence of the 
limit Coo °f the sequence of real constants Cj, connected with <f>j, 
either finite or infinite. 


Thus 


Coo = lim C 3 

j^oo 


(42) 


must exist, either finite or infinite, provided the phase-locking takes place. 

On the contrary, the phase-locking does not arise and the phase function 
reveals apparently ‘irregular’ behavior if and only if the sequence Cj has no limit. 


Remark: 


Allowing constants Cj and their limit to assume infinite values, the 
real projective line RP 1 isomorphic to the circumference has to be 
adopted as the space of their (C-constants) originating. This can be 
realized by means of identifying each Cj with the pair (Cj, 1) and 
the interpreting the latter as homogeneous coordinates on RP 1 . The 
limit of the sequence Cj has also to be understood as a point in RP 1 . 

Figures [2H1 EH display the two examples of C-sequences for the phase-locking 
(Fig. HH) and ‘irregular’ (Tig. I3T1) phase evolutions, respectively. The distinction 
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Figure 21: Here the sequence Cj not tending to an apparent limit indicates the ab¬ 
sence of the phase-locking (‘irregular’ phase behavior). The bias is the same periodic 
rectangular pulse sequence as above besides the specific value of the DC contribution 
Me = 1-40. 

arises because to the slightly different values of the constant bias parameter: 

Me = 1-46 for Fig. [2H1 and t dc = 1.4 for Fig. ETJ cf. the left-lower panel of Fig. [TUI 
The other parameters of the bias function / are the same as for the pulse shown 
in the inset in Fig. [U| 

In order to detect the phase-locking without a plain phase simulation, one 
may, in spite of determination of the constants Cj from the functions 4>j(t) by 
means of Eq. (EH, make use of the following recurrence relation which can be 
derived from Eqs. EH, EH, EH and Cj definition: 

_ Cj (e _ ^ p (°>( T ) cos 2 <^ ( 0 ) (T) — Q ( 0 ) (T)e^ p (°)( T ) sin § 0 (o)(T)^ — e^ p (°h T ) sin 50 (o)(T) 

Cj (e“ 5 p (o)( r ) sin§0 (o) (T) + Q { 0 )(T)e^ p (°) (T) cos |0( O )(T)) + e5 p (°) (T) cos ±0( O )(T) 

(43) 

(It is worth reminding that the subscript ‘( 0 )’ hints at the specific initial condi¬ 
tions the functions 0(o), P( o), Q{ o) obey which read 0(o)(O) = 0 , -P(o)( 0 ) = 0 , Q(o)( 0 ) = 

0. It emphasizes the distinction of the ‘ground’ solution and its derivates from the 
particular element 0 o of the sequence 0 ^; we shall omit this notation complication 
whenever possible.) The relation above is of a notable importance allowing one 
to obtain, finally, the exhaustive description of the long-term phase behavior on 
the time scales exceeding T by means of pure algebraic manipulations. 

To that end, one has to mention that the map Cj —> Cj + i implied by Eq. (1-iTfli 
is fraction linear. It is well known that any such transformation is equivalent to 
a linear automorphism on a projective space. To make use of much simplification 
following from such a linear reduction of the problem, let us consider the vector 
space of 2-element columns C = (using here and below boldface characters 
for notation of matrix-valued quantities) with the ratio of the elements equal 
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to C, a/b = C. We consider the columns differing by a non-zero multiplier as 
equivalent, [^] ~ [^?] , A ^ 0. Let us also introduce the 2x2 matrix 

^ = ( e~5 p o( T ) COS §0o(T) - Qo(r)e3 p o( T ) sin i^j(T) -e5 p °( T ) sin ^(T) \ 

<l> " \ e-5 p o( p ) sin I0 O (T) + Q 0 (T)e5 p °( T ) cos ±0 O (T) +e^ p °( T) cos I</» 0 (T) J 

(44) 

which is, as a straightforward check shows, unimodular, det = 1. Then it is 
also straightforward to show that the transformation 


Cj+1 = 4>C, (45) 

is equivalent to Eq. (@3D- The map C j i—>• C,-+i corresponding to Eq. P^jl is just 
the linear automorphism of the real projective line RP 1 mentioned above. 

Since the transformation associated with matrix $ does not depends on j, 
starting from j = 0, after j — 1 iterations one comes to the equation 

Cj = &C 0 . (46) 

It is convenient to chose C 0 = [ j 0 ]. The very ‘starting’ C-constant C 0 encodes 
the initial value of the phase function <p{t). Giving 0(0), it can be calculated by 
means of Eq. (I2U1) . 

Eq. (14611 is, essentially, the desirable algebraic relationship describing the long¬ 
term phase evolution. The time variable t is here encoded in the integer variable 
j playing role of a discrete time on the scale T and amounting, numerically, to 
the integer part of t/T. 

We shall derive below the expanded version of (USD where the power of the 
operator $ is given in explicit form. 


The ‘irregular’ phase evolution 

We consider, first, the case of ‘irregular’ phase evolution when the discriminant 
A, Eq. (1401) . is negative. In this case we proceed with the problem of calculation 
of the accumulated phase variation over a long time interval (the single period 
averaging does not yield a meaningful result here). In accordance with (jJJ) , it 
immediately yields us the value of the average voltage across the junction, i.e. 
the physically measurable quantity. In particular, we shall see that, in spite 
of apparently ‘irregular’ phase behavior, the phase evolution manifests actually 
no signs of chaos. Moreover, as a matter of fact, the phase reveals, in a sense, 
oscillation of a definite frequency and the apparent irregularity of its evolution is 
manifested because this frequency is in general incommensurable with the bias 
guiding one (quasiperiodic behavior 0). Asa particular consequence, the average 
voltage converges to a quite definite value which can be calculated in a general 
case. 
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Elaborating the relationships referred to above, let us calculate, at first, the 
eigenvalues of the matrix In view of its role in the phase evolution, one should 
not be surprised that the discriminant of its characteristic equation 

A 2 + 1 - 2A (cosh ± P 0 {T) cos ±0 O (T) - ^ Po(r) Q 0 (T) sin ±0 O (T)) = 0 (47) 

coincides, up to the positive factor Ae~ Po{ - T \ with A = A[/], see dim . Thus, it 
is precisely the case of the ‘irregular’ phase evolution when the matrix flU has 
a pair of complex (and complex conjugated) roots A±. Furthermore, since their 
product is the unity, they equal exp ±ia for some real a. More concretely, one 
gets 

A± = e ±ia = (48) 

i e iPo(T) (_g o(T ) sin 1 0 O (T) + (1 + e ~ p ^) cos ± <h{T) 

± iv 73 ^) • (49) 


These eigenvalues are connected with the following (complex valued) eigenvectors: 


V± = 


-2 sin ±0 O CO 

Qo(T) sin i MT) + (1 - e- p » <T) ) cos i fo(T) ± iv^ 


(50) 


In other words, 

<£V ± = e ±ia V ± . 

Further, the following identical decomposition of an arbitrary real-valued 2- 
elcment column A in terms of the columns \± takes place: 


A = 



K+{ A) V+ + AT (A) V, = 2ft(K+{A) V+), 


(51) 


where we make use of the notation 


x 


A±(A) = ±i 4sin±0 o (Ov-A 

A (Qo(0 sin ±0 o (O + (1 - e _Po(T,) ) cos ±0 o (O T iv 7 ^ 

+2 B sin ±0 O (T) 


(52) 


The j-fold (j = 0,1,2,...) application of the linear matrix operator $ to the 
vector A expanded in accordance with (EH) leads to the equation 


A = e lja K+ (A) V+ + e~ ija K_ (A) V_ . (53) 


This is in fact the desirable explicit representation of the operator power 
convenient for our purposes. 
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Applying the decomposition (E2) to Eq. (1551) . the resulting equation can then 
be resolved with respect to Cj representing it as a fraction-linear function of Co, 
Its formula will be given later on (see Eq. (HI) while here we write down instead 
the explicit form of the equation 


e i 4>(jT) = i<fe(0) = gi^nhO) 1 + Q-^(0)(°) = 1 lC j 

1 + CjJ~( o)(0) 1 + iCj 


(54) 


following from the definition of Cj (and making use of the initial values </>(o)(0) = 
0, JF( 0 )(0) = i, following, in turn, from definitions of </>( 0 ), -A(o))- It reads 


e W 


e h Q C(+) + 

----, j — 0,1,2,..., 


where all the coefficients 


(55) 


f/(+) = riR + v^A - i (ru + CqV^a) , 

U^ = -n R + V-A + i - CoV-A^ , 
n R = C 0 (l - e' p (°) (T) ) cosi0 (O )(T) + (C 0 Q(o)(T) + 2) sini0 (o) (T), 
ru = (1 - e _p (°) (r) + 2C 0 Q( 0 )(T)) cos §</> (0 )(T) 

+ (2 C 0 e~ p ^ + Q m (T)) sin |0 (O) (T) (56) 


do not depend on j and n R ,nj are real. (Notice that the only ‘complex valued 
ingredient’ in is the factor i.) Since j represent here the time ‘normalized and 
discretizied on the scale T' (i.e. the integer part of t/T, in fact), Eq. (1551) manifest, 
essentially, the specific ‘hidden’ periodicity of the phase function which has been 
occasionally referred to above, the corresponding period amounting to 2nTa~ 1 . 
Generally speaking, the latter quantity is incommensurable with the bias period 
T. As a consequence, the discrete parameter j never assumes sequential values 
differed by 2nTa -1 or a quantity aliquot to it. It is this circumstance which does 
not allows one to distinguish the oscillations which would be described by Eq. 
1551' if j were a continuous variable. 


Average rate of the phase growth 

Now we are in position to apply the relationships derived above for the determi¬ 
nation of the average voltage I4 V across junction in the case of ‘irregular’ phase 
behavior (phase-locking absence). In view of Eqs. (J2J), (J5J), 

¥A« = I«jT)-#0)] (57) 

J 
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and, therefore. 


exp 


2br T 

% 


V U) 

av 


{e^ jT ) e 



1 + iC 0 1 7 f e ija U^ + 

1 - iCo J \ e -ijajj(+) + e ijaf/(-) J 


(58) 


Thus the averaging over a lapse consisting of unboundedly increasing number 
of bias periods reduces to the computation of the limit V= lim^oo for 
V// satisfying (jHHl) . It yields a definite result if and only if the latter exists. Let 
us consider how it can be computed. 

As to the first multiplier involved in Eq. (USD, its nonzero base does not depend 
on j and the limit always exists (obviously, it equals the unity). On the contrary, 
the second multiplier is not a universal constant. It depends on the ratio of the 
modules of the complex coefficients U and U^~\ 

Clarifying the latter point, let us calculate the difference £ = | UC) | 2 — | JJC) | . 
A straightforward algebra yields 


£ = 8a/— Ae P(0) £o, where 
£ 0 = C„ (e p wQ { o) cos ±0 (O ) + sin ±0 (O )) 

+C 0 ((e P(0) - l) cos |0(o) + e P(0) Q( 0 ) sin ±0 (O )) + e p w sin |0 (O ) (59) 

(the argument (T) is here omitted, cf. (15011 ). The sign of the difference coincides 
therefore with the sign of the factor £ 0 , Apart from the constant parameters 
determining <j>(o), P(o),Q(o), S 0 also depends on C 0 which encodes the starting 
(initial) value of the phase and may assume, in principle, arbitrary real value 
including the limiting case C —> oo. However, £o is itself a second order polyno¬ 
mial in C. Calculating its discriminant, one finds it to be connected, again, with 
A [/], equaling to 

e 2Po(T) A, (60) 

and being in our case negative just in view of the ‘irregular’ phase behavior 
condition A < 0. Therefore, as a function of Co, £o has no real roots and 
is either everywhere positive or everywhere negative irrespectively of the initial 
phase (connected with the specific value of Co). Furthermore, since 

£ 0 |c 0 =o = e Po(T) sin |0 O (T), (61) 


the sign of £o (and £) coincides with the sign of sin which is determined, 

ultimately, by the function / and may not vanish by virtue of definition of A and 
the ‘no convergence’ condition A < 0, see flUl). 

Thus we have shown that 


|f/(+)| > | £/(-) | iff sin I0o (T) > 0 
and 

\U(+)\ < \uC)\ iff sin 100(T) < 0. 
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One may employ any of the two decompositions shown below of the second 
multiplier in the last line of 


g i joLjj{+) _ g-ij l v 


e-ijaJjM - e ijay(-) 

f/(+)' 1 




2ia 


—2i a 


f/(+) 

f/(-) 


1 - e~ 2ija 

1 — e 2 ««(*/(-)/*/(+)) 
1 — e 2 b« ([/(+)/[/(-)) 
[/(-) 7 I 1 - e -2iia(f/(+)/[/(-)) 


(62) 


Let us arrange to apply the upper representation if |tA + l| > and the lower 

one in the opposite case. Then the limit of the last multiplier in the corresponding 
line of (1531) exists and are equal to the unity since its base to be raised to the 
power 1/j is a finite ^-independent gap distant from zero. The existence of limits 
of the second multiplier is evident, it also equals the unity. Thus the limit of the 
first (constant) factor yields the limit value for the whole product. 

Having thus calculated the limit of una as j —> oo, Eq. (EED leads to the 
following conclusion: 


exp 


2br T 
<h 0 




e 2ia for sin |0o(T) > 0 
e~ 2ia for sin I0o(T) < 0 


(63) 


where a is defined in Eq. (1481) . This simple result yields the explicit representation 
of the average voltage across junction in the case of ‘irregular’ phase evolution: 


$0 f ol + kit if sin |0 O (T) > 0 
7 tT X \ —a + k'K if sin|0 o (T) < 0 


(64) 


for some integer k (we shall discuss the method of its determination later on). 

More exactly, the formula above describes, for a single k, a single branch 
binding the two neighboring Shapiro steps, i.e. ‘horizontal’ constant voltage seg¬ 
ments on the junction I-V curve whose orders differ by the unity. Such branches 
are sometimes called ‘resistive portions of I-V curve’ |Sj although the specific 
dependence of the average voltage on tdc, as it is described by Eq. dOil), may 
considerably deviate from the simple Ohm’s law proportionality. In particular, 
near the edges of these ‘resistive portions’ the differential resistance dV^ 00 ^Vdtdc 
diverges. 


Solution of the spectral problem in the case A > 0 

After a minor modification, the majority of the above equations is also applicable 
in the case of phase-locking. Here one has to assume 

A > 0. (65) 
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Then Eqs. (14811491) hold true while the adapted versions of Eqs. (©-© read 


A± = |e.2 Po(T) Q 0 (T ) sin ±<f> 0 (T) + (1 + e Po{T) ) cos ±<po(T) 

± \/A) , (66) 

V = r —2 sin \(j) 0 {T) 1 , fi7 * 

V± [ Q 0 (T) sin I0 O (T) + (1 - e“ p °( T )) cos §0 O (C ± J ' '• ] 

These are the eigenvalues and the eigenvectors of the matrix $ (14411 . and the 
following equations take place 


$V ± = A±V±, 


where in this case A±,V-t- are real. It is worth noting that the unimodularity 
constraint 

A + ■ A_ = 1 (68) 

also holds true. Then it follows from Eq. (14611 

- 4 sin i0o (T)VACj = X\K + V+ + A j _K_ V_, (69) 

where 

K± = K±(C 0 ) =C 0 L^+ 2 sin ±MT), (70) 

L± = Qo(T)sini0o(T) + (l-e- Po(T) )cosi0o(T)±v / A, (71) 


and, after some algebra, one gets 


Cj = —2 sin 300 CO- 


KK + (C 0 ) + \ j K-(C 0 ) 


XiK + (C 0 )L + + X J _K_(Co)L_ 

-2sini0 o (T) x 

(L_C 0 + 2 sin i0 o (r))A J + - (L + C 0 + 2 sin i0 o (r))Aj 
(L-C 0 + 2sini0 o (T))L + Ai - (L + C 0 + 2 sin ^ 0 {T))L_X j _' 


(72) 


Here all the dependence on j (proportional, up to normalization and dis¬ 
cretization, to the evolution time t) is isolated in the factors A0. More precisely, 
the j dependent terms combine to the powers (A+ZA.)- 7 or (A_/A^)- 7 but, in view 
of Eq. m, these coincide with A±, respectively. 


Calculation of phase function 

Basing on (ED, the explicit equation completely determining the phase at any 
moment of time in terms of the ‘ground’ phase function 0(o) specified on the 
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interval [0, T) and the functions P(o), Q( o) calculated from 0( O ) and defined on the 
same interval is easily derived: 


e i 4>{t'+jT) _ gi0(o)h') x 

(■ M + L + Xj - M_L_XL) - 2sini0( O) (r)(M + A J + - 

(M+L +- M_L_ XL) - 2 sin §0 (o) (T) (M+ - M_ Ai)P (0) (t') ’ 

= x 

(£+ - 2sm^, cl) (r)J^(i 7 ))Af + A J + - (L_ - 2 3 in^„, ) (r)7i^(i ; ))M-A J l 
(i+ - 2 s ini0 (o) (r)^ (o) (i'))M + Ai - (L. - 2 s ini^ 0 ,(T)J' (0) (i'))M_Ai A 

where 


M ± = M ± (C'o) — h T C 0 + 2 sin |0( O )(T), (75) 

t' G [0, T) is the excess of t over the nearest lower jT, t' = t — T[[t/T]], for integer 

J = M. 

Remark: 

Eqs. (d-d also apply to the case of ‘irregular’ phase evolution and 
exactly in this form. The only difference is in the definitions of A±, L± 
©,CID- In ‘irregular’ case one has to replace the real term y/A by 
the pure imaginary i\/—A. Then {A ±, L±, M±) = {ALzp, M T } that, 
in particular, ensure the phase function defined by JED to be real. 

Eq. (1731) determines the phase at any moment of time (up to a constant 
summand aliquot to 2ir which will be computed later on) through the phase 
initial value 0(0), the solution 0 O of Eq- ® vanishing at t — 0 and the complex 
valued function P 0 (f) both computed on the finite interval [0, T]. The exhaustive 
description of the process of the asymptotic establishing of the steady phase¬ 
locking state is its particular byproduct. It allows one to explicitly describe the 
‘transient processes’ representing the distinction of arbitrary given phase function 
from the ‘refined’ steady one to which the former asymptotically converges. 

To be more specific, let us now consider in brief some straightforward conse¬ 
quences of Eqs. jnikfzsD- 

At first, it should be noted, that the dependence of r.h.s. of JED on j (the 
time normalized to the scale T and then ‘discretizied’) disappears in the case of 
satisfaction of any of the two equations 

M+ = 0 & C 0 = -2 sin I0 O (E)/L_ or (76) 

M_ = O^C 0 = -2 sin I0o (T)/L+. (77) 

(It is worth noting that, requiring Cq to be real, each of them implies the condition 
A > 0. Thus, neither of Eqs. ClSD, jzzd can be fulfilled in the ‘irregular’ case.) 
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Then e 1 ^ is a periodic function of t with the period T. In particular, it coincides 
with the own asymptotic limiting form. Such phase evolution can be named a 
steady one. Notice that the ‘steadiness’ does not means the periodicity of <f>(t) 
but allows, apart from periodic part, the uniformly growing contribution 2tt kt/T 
for an integer k. 

It can be shown by a direct check that C 0 , satisfying m or (EH), also verifies 
Eq. (1551) . Hence, (TTUj) and (EZI) determine the C-constants which correspond 
to initial data characteristic of the ‘refined’ (steady) phase-locking states. Two 
choices ± correspond to two states, stable and unstable. Such C-constants were 
referred to above as C^. Now let us arrange to reserve this notation for the 
C-constant implying the stable steady phase evolution. 

In order to determine which of these two states is stable, let us consider a 
generic case when both M + 7^ 0 7^ M_. It follows from (16611 and (16811 that the 
condition A > 0 implies either |A + | > 1 and |A_| < 1 or |A+| < 1 and |A_| > 1 
(and both A± are either positive or negative). 

Let us arrange about the following interpretation of the subscripts 
‘max ’ and ‘min 


if |A+| > 1 & | A_ | < 1 
if | A_| > 1 & | A+| < 1 


(78) 


then ‘max ’ means ‘— ‘min ’ means ‘+ ’. 


Then one may introduce the equation 


—Qo(T) sin §</> 0 (T) + (1 + e Po(T) ) cos ±0 O (T) | - \/A 
~Qo(T) sin I0o(T) + (1 + e- p °( T )) cosi0 o (E)| + v^A 

(79) 

defining by them the positive constant x. Using it, Eq. El can be recast to the 
following form: 


2 x _ /K min 


e i Ht'+jT) 
where Z(t) 


JM?) Lmax ~ 2sini0 o (r)Jo(f) 1 - Z(t')e 2 ^ 

' L max - 2 sin \<j> 0 {T)F 0 {t') ' 1 - Z(t')e~^ ’ 1 J 

Lmin 2 sin i0 o (T)jEo(t) M min 
L max - 2 sin i0o(T)J^o(f) M max 


Sending here j —► 00, one hnds that the exponent e 1 ^ for a ‘generic’ solution 11 
4>(t) of Eq. © approaches to the periodic function edefined for t e [0, T] by 
the equation 


e i0oo(i) _ gi^oOO . ^ max 2 sill 10o(r)jF o (f) 

Lmax 2smi0o(r)J r o (()' 


(81) 


11 for which Z(t) ^ 0 
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Notice that the mentioned periodicity is not manifest from this formula, it is 
rather a specific consequence of the properties of its ingredients constituting the 
essence of the definitions of the latter. 

The very function of the ‘refined’ phase-locking steady phase evolution, (j)^ (£), 
advances at each time step of duration T by 2irk for some fixed (t-independent) 
integer k. The formula m is nothing else but Eq. ( 12 % 1 ) with 


C 


C^ = 


2sinf</>o(T) 


L 


max 


(82) 


This answers the question which of the two equations M + = 0 or M_ = 0 
describes the stable steady evolution (asymptotic phase-locking state): it is de¬ 
scribed by the equation 

M min = 0. (83) 


Phase-locking criterium 

The detailed picture of the convergence of a generic phase function fi(t) to the 
asymptotic limit 0^ (t) can be inferred from the equation (1KU1) recast as follows: 


e i m'+jT)-^(t'+jT)) 


1 - Z(t')e~ 2 ^ 

1 - Z(f)e- 2 ^’ 


t' e [o ,T),j 


0 , 1 , 2 ,... 


(84) 


which describes what can be called a ‘nonlinear exponential’ convergence. If 
j = [[t/T]] is sufficiently large to ensure the satisfaction of the condition 


max \Z\ < e 2 ^- 7 , (85) 

[0,T] 


the simple estimate follows 

|0(f / +jT)-0 oo (f / +jT)| < log(l + \Z\e^ 2 ^) —log(l — \Z\e~ 2 *i) < 2max|Z|e~ 2 ^ 

which describes just the exponential convergence. The factor max|Z| and the 
exponent coefficient x both determine how many periods (enumerated by the in¬ 
teger j ) have to elapse until the phase function approaches the steady asymptotic 
state with the prescribed accuracy. 

At the same time, for comparatively small time t, the ‘transients’ 0(t) — 0oo (t) 
may have no relation to the exponent. Eq. (USD allows to estimate the duration 
of this ‘near-zone’ lapse. 

As to the second possibility 

ri - 2sini0 o (T) 

^ — ^tX] - j 

^min 
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(the solution of the equation M max = 0), it also yields, formally, the steady phase¬ 
locking state described by the phase function 0 M (t), which can be computed with 
the help of equation 


gi 0x(i) _ gi0o(i) . Lmin 2 sill |(/> 0 (T) JFq (t) 

L min - 2 sin \(f 0 (T)F 0 (t) ’ 


(87) 


and is distinct from 4>oo{t). This is unstable solution ‘repelling’ any neighboring 
one. Indeed, having rewritten Eq. PI) as follows 


gi <t>(t+jT) 

where Z(t) 


„*.<« . 1 - 

1 - Z(t)e 2 "-> ’ 

Lmax 2sini0 o (r)J'o(f) M max 

L min - 2 sin \<j>o{T)F Q (t) M min 


( 88 ) 


one sees that for arbitrary small but non-zero M max (which is non-zero if <f> ^ 
mod 27 t) the function e 1 '^7+J T ) escapes exponentially off the function e 1< Axb+j' r ) 
which, itself, can be revealed only in the case of the exact vanishing of M max . For 
any other phase function ultimately approaches, as j increases, to e 1 ^ 00 . 

The computations above proves the following statement which has been men¬ 
tioned above: 


Theorem 1 The condition A > 0 is necessary for the phase-locking to be ob¬ 
served whereas the strict inequality A > 0 is sufficient. □ 

Remarks: 


• There is a formulation of the criterion above operating with not specific but 
arbitrary solution <f of Eq. © (instead of </>(o)) and not requiring for the 
functional T to obey the specific initial condition CU) 12 . In general form, 
it reads 


\D\ > 

where D = 
A(f = 


1, 

3[e~ W(;T(t 0 + T/2) -ggjj - T/2))] 

2 V^(to + T/2)]S[.F(to - T/2)\ 
4>(to + T/2) — f>(to — T/2). 


D does not depend on f 0 . 


(89) 

(90) 

(91) 


• In terms of D, the eigenvalues A± are represented as follows 

A ± = D± VD 2 - 1. 

• It is shown below that for A = 0 the phase-locking is also observed but its 
properties reveal some distinction from ones of the case A > 0 (the ‘weak’ 
phase-locking against the ‘generic’ one). 

12 Nevertheless, T(0) may not vanish for any representation. 
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Weak phase-locking 


For the sake of completeness, let us consider here the specialities of the situation 
intermediate between the exponentially stable ‘generic phase-locking’ taking place 
if A > 0 and ‘irregular’ phase evolutions for which A < 0. ft is distinguished 
by the merging of the two eigenvalues of the matrix (jTflh Zero discriminant 
condition and Eq. © imply 

cosh \Pq{T) cos (T) — fe5 Po(T) Q 0 {T) sin \(j>o{T) = either 1 or — 1 

and then it follows from (JHHI) that the two-fold eigenvalue of $ which now is 
represented as follows 

A = i e 2 p °( T ) ((1 + e~ p °( T )) cos ±<f) 0 (T) — Qo(T) sin §</> 0 (T)) (92) 

has also to be equal to either 1 or - 1 . 

In the case A = 0 Eq. (EH) does not apply. Instead, one may employ the 
following explicit representation of the power of the 2 x 2 matrix 

* = (“$) (93) 
whose elements obey the constraint 


46c + (a — d ) 2 = 0 


(94) 


(just meaning that the $ eigenvalues coincide) times the 2 -element column with 
arbitrary elements A, B : 




A 

B 


(i {a + d)Y 


B 

2 c 


a — d 
2 c 





+ 


a -\- d 



(95) 


(It can be established, for example, by means of the mathematical induction.) 
Since the eigenvalue does not vanish, a + d ^ 0. Let us assume also, for a while, 
that c ^ 0 / ft which imply, in view of (1941) . a — d ^ 0. Then the following 
analogue of Eq. (1721) arises: 


_ (a + d)Co + ((a — d)Co + 2 b)j 
J (a + d) + (2cC 0 — (a — d))j 

It implies the existence of the limit 

lim C _ ( a ~ d)C 0 + 26 _ ( a ~ d)Co + 26 
j^o o 1 2cC 0 — (a — d) AbcC 0 — 26(a — d) 
2 ^ (o — d)Co + 26 —26 

(a — d)((a — d)Co + 26) a — d 
a — d 
2 c 


(96) 


(97) 
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which proves independent of Cq. It also follows from © that Cj does not depend 
on j (all the elements of the sequence coincide) if and only if 


C 0 = 


Coo 


—2b 
a — d 


a — d 
2 c 


(98) 


Thus, as opposed to the ‘generic phase-locking’ taking place for A > 0, there 
is only a single initial phase which yields the ‘refined steady’ evolution. All the 
other phase evolutions converge, with the course of time t, to the latter which 
plays therefore the role of attractor. 

Interestingly enough, the same phase function plays simultaneously the role of 
the repcllcr and all the distinct phase functions approaching it are simultaneously 
‘moving away’. There is not contradiction here since the phase functions live 
in fact on the closed circumference and ‘the one side attracting’ to a point is 
simultaneously ‘the another side repelling’ from the same point. 

To illustrate the simultaneous attraction/repclling property of a steady phase 
function, let us compute the leading j-dependent contribution to Cj. The result 
is 

Cj = C. + + o(r 2 ) (99) 

The terms shown do not depend on the initial phase (encoded in Co) which af¬ 
fects only the higher order contributions. This means, in particular, that all the 
phase functions approach their common steady limit ‘from one side’ (the leading 
contribution to the deviation from the limit is common for all of them). Then, 
obviously, the closer the initial phase is to the one corresponding to the refined 
steady state from this side, the less time is necessary for the reaching, in appro¬ 
priate sense, the steady limit. On the other hand, the ‘very long’ phase evolution 
until it reaches some fixed vicinity of the ultimate steady state arises when the 
initial state is ‘very close’ to the steady phase function from the opposite, ‘wrong’ 
side. This behavior differs from the ‘generic’ phase-locking where the choices of 
the initial phases closed to the phase of the stable steady evolution always lead 
to the quick ‘monotonous’ convergence to the limiting function irrespectively to 
the initial relative angular direction. 

This specialty can be inferred more rigorously from the consideration of the 
analogue to Eq. ( 171 ) which follows from © and now reads 

gitfi'+Tj) = e i 4>( 0) (t') x ( 10 0) 

2HG + (l + CoT^d)) + (G_ + 2HC 0 ){2H - 

2HG + (l + CoT(g){t')) + (G_ + 2HC 0 ){2H - ’ 

where 

G± = (1 ± e _p w (T) ) cos |0(o) (T) T Q(o)(T) sin ±0 (O )(T), 

H = e _p w (T) sin §0 (o ,)(T) + Q(o)(T) cos i <f>( 0 )(T) (101) 
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(note that, in particular, G + = 2e~^ p ^ 0) \ equals either 2 e“ 2 P (o) Q r —2e _ 2 P (°> and 
does not vanish). It determines the phase function for arbitrary t expressing it 
through the functions 4>(o), P{ 0 ), Q(o) specihed on the segment [0, T], The most 
substantial difference with Eqs. dSSD and (IZSD is the non-exponential dependence 
on j. Now the convergence to the asymptotic phase function </>“(i), defined 
mod 27r for any t by the equation 


gCW = gi^(o) (t 1 ) 2H G-Fwit) 
2H-G-F i0) (t') 


t’=t-T[[t/T}\ 


is linear in j 1 oc t 1 . 

Eq. (11001) can also be rewritten as follows 




where Z(t') 


~Z{F)SC + 1 + 2 HG+'jSC 
Z(t’)SC + l + 2HG+ 1 jSC' 

_ G_ 

2 H-G.T(t'Y ' ' 0+ 2 H' 


( 102 ) 


(103) 


For bounded j, the arbitrary constant Co can always be chosen making SC so 
small that the unity is the dominating contribution in the both numerator and 
denominator in m- It makes their ratio to be as close to the unity as one 
desires. If further j unboundedly increases, depending on the sign of SC, the 
two distinct situations can occur. Namely, if SC is of the same sign as HG+ 1 , 
increasing j, the ratio monotonously’ tend to the unity, the less SC, the 

faster the limit is reached. This is the case of the ‘true’ situating of Cq with 
respect to initial phase C*, of the steady phase function. However, if SC is 
still small but has the sign opposite to the sign of HG+ 1 , increasing j, the j- 
dependent contribution is subtracted from the leading terms of the numerator 
and denominator (here the unity) and the absolute value of their sum decreases 
reaching the minimum for j = [[|(//hC)” 1 G + |]]. Simultaneously, the ratio (II Uhl) 
varies in some way and approaches the initial closeness to the unity only for 
j ~ 2\(HSC)~ 1 G + \ or greater. The further increasing of j already leads to the 
‘monotonous’ converging to the unity with the rate ~ j -1 as above. However, the 
‘convergence time’ estimated as 2T\(HSC)~ 1 G + \ increases as the deviation SC 
of the initial phase from the phase of the steady function, having ‘wrong’ sign, 
becomes smaller. 

Finally, to abandon the temporal assumptions made above, it has to be noted 
that the cases c = 0 or b = 0 implying a = d which are not covered by the formulae 
above, reveals qualitatively similar relationships. Their analysis is carried out in 
a similar way, provided the following simple representations of 


a 
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0 

a 


= a 3 


1 M 

jc/a 1 J ’ 
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are utilized. They also lead to the convergence of the order ~ j~ x . 

Resuming, in the case A = 0 not covered by the criterium formulated in 
the form of the Theorem |U the phase-locking understood as the asymptotic 
convergence of all phase functions to some steady one reproducing the own form 
mod27r on each segment of t variation of length T takes place as well. However, it 
proves not exponential in time. Rather, the steady limit is approached as ~ t" 1 . 


Winding index 

Let us now derive yet another important property of phase-locking solutions of (l9l) 
following from Eq. (HU). It enables one to calculate the integer k entering, in par¬ 
ticular, equations d3U><inu- To that end, calculating the logarithmic derivative 
of Eq. (EBP and applying (EU, one gets the equation 


i d 0 = id 0 (o) + 


CdL 


(o) 


C d T, 


(o) 


1 + CjF(o) l + CT(p) 

d t e^(°) d t 
— id0(o) + i C - —-b iC 


= i d 0(o) + 2\CTt 


1 + CT (o) 

e i ^(°) 


1 + CT(Q) 

Integrating it on the interval [0, T], one gets 

0(T)-0(O) = 0 ( o ) (T) + 2CK 


1 + CT (o) 

d t 


r 

e i0(°) 

L 

_ 1 + CT (o) _ 


d t. 
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Up to this point, the manipulations above do not go beyond Eq. m, definitions 
and identities. The situation drastically changes if we substitute here in place 
of arbitrary C the constant C^ (1821) which, in the case A > 0 here assumed, 
converts an arbitrary solution 0(f) of Eq. © to the stable refined phase-locking 
phase function 0^ (f). Since the latter advances on each time step [f,f + T] by 
the strictly fixed increment 0oo (f + T) — 0oo(^) = 27Tfc, one gets 

Theorem 2 Let 0(o) be the solution of 0) obeying the initial condition 0(o)(O) = 
0 and determining the complex valued function T which, in turn, satisfies Eq. 
m> with initial condition JF( 0 )(0) = i. Let also A[/] > 0, where A[/] is defined 
by Eq. !&$) . Then 13 
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-sin i0(o)(T) 


St 


7T 


gU (o) (0 

2 sin 10( O ) (T)T ( 0) (f) _ 


d t 
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is the integer equal to the total sum (taking into account the orientation sign) of 
the number of full revolutions (the winding number) which any phase function, 

13 There is a misprint in the corresponding equation in [JJ. 
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except of 0c*], will carry out on the time steps of duration T in its asymptotic 
steady phase-locking state. 

Remarks: 

• The formula mm provides us with the substantiation of the statement 
made above which concerns the constancy of the order of phase-locking 
when the bias parameters varies. Indeed, in view of (EH) the order k is 
continuous with respect to the variables parametrizying the bias function 
/. In particular, k continuously depends on tdc as far as A retains positive. 
Hence it assumes a constant value on the i ( \ c intervals over which A graph 
shown in Fig. El is situated above the horizontal coordinate axes (more 
generally, as the bias function / varies within the connected component of 
the phase-locking area). 

• A more general representation of k which operates with not special but 
arbitrary solution of Eq. (El) can be derived from Eq. (HUD , cf. the re¬ 
mark following Theorem 1. However, as opposed to the case of criterion 
based on the ‘universal’ expression EH, here the two inequivalent formulae 
distinct by the opposite roles of the left and right boundary points arise. 
Their difference is a complicated non-linear (and apparently non-trivial) 
functional which vanishes on all solutions of El), provided the bias function 
corresponds to the phase-locking phase evolution. 


Conclusion 

The equation El) and the properties of its solutions seem to be of a considerable 
interest in view of several reasons. First of all, this is, of course, their sound 
physical relevance following from the extensive applications in the applied theory 
of electric activity of Josephson junctions which employs Eq. © as the base of 
the efficient model of the junction phase dynamics m in the important case of 
negligible role of junction capacitance (overdamped junctions) On the other 

hand, Eq. © is of evident interest in its own rights from a pure mathematical 
point of view. It provides us a remarkable example of apparently supreme simple 
non-linear ODE which prove associated with a linear problem 14 and, in view 
of such a link, allows a deep exploration by analytic methods. At the same 
time, in spite of its apparent simplicity, it is definitely far of being regarded 
as a mathematically trivial entity. It suffices to say that the properties of Eq. 
El) are still not completely understood even for sinusoidal bias function / = 
a + bsm(cut + to), the case of a primary interest from viewpoint of applications. 

The relationships considered above does not exhaust the collection of rigorous 
ones which © allows to establish by means of elementary technique. However 

14 This point mentioned also in footnote 171 InagelT^ll is beyond the scope of present discussion. 
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they are distinguished by the advantage of universality being valid with fairly 
weak limitations on the class of allowable bias functions. 

The introduction of the functional A[/] (1401) is the central point of the ap¬ 
proach. Eq. dZSD is noteworthy as the explicit representation of the phase function 
for arbitrary t through a single solution of Eq. © computed on the finite interval 
[0, T\. It is this equation which allows to establish the convergence, in the case 
A[/] > 0, of any phase function (except of see Eq. (ED) to an asymptotic 
limit and to show that this limit coincides mod 2n with (f)^ defined by Eq. ED- 
Thus Eq. (ESP yields the rigorous model of the phase-locking property allowing 
one to compute any of its quantitative characteristic of interest. 

In combination with Eq. ED describing the long term phase evolution in 
the opposite case A[/] < 0, the above relationships lead to the criterion of the 
asymptotic property of the phase-locking (Theorem 1) which can be reformulated 
to operate with arbitrary single solution of Eq. © on a finite segment of the length 
T and its derivates. 

Another important result is the formula m determining the integer winding 
number k (phase-locking order) through the same easily computable data. In 
physical terms, the latter integer quantity is directly connected to the average 
voltage applied across a junction in the phase-locking state which is the supported 
constant and proves independent of the slow variations, up to a certain extent, 
of the parameters, provided the period T is kept unchanged. This effect lies in 
the core of the modern DC voltage standards |5j. 

In view of the above remarks, the properties of Eq. © and its solutions is a 
fruitful area of the mathematical study which is worth of a further development. 
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